Major questions:
Firstly, we’re going to install and load our packages. We will use the pacman Package Management Tool which allows us to install and load subsequent packages in a condensed and efficient way.
# Install and load required packages
pacman::p_load(dplyr, ggmap, ggspatial, ggplot2, janitor, lme4,
readxl, sf, tidyverse)Secondly, we’ll read in our data and explore it.
Let’s view the data.
Now, we’ll wrangle some of our data to get it ready for analyses and plotting.
The weather data were acquired from the Australian Bureau of Meteorology weather station 70351.
# Subset to actual instances of bait take by a vertebrate
data <- data %>%
subset(species != "Lizard") %>%
mutate(fox_visited = ifelse(fox_visit_1 != 0, "Yes",
ifelse(fox_visit_2 !=0, "Yes",
ifelse(fox_visit_3 !=0, "Yes",
ifelse(fox_visit_4 !=0, "Yes", "No"))))) %>%
mutate(fox_visits = ifelse(fox_visit_4 != 0, 4,
ifelse(fox_visit_3 !=0, 3,
ifelse(fox_visit_2 !=0, 2,
ifelse(fox_visit_1 !=0, 1, 0))))) %>%
mutate(raven_visited = ifelse(raven_visit_1 != 0, "Yes",
ifelse(raven_visit_2 !=0, "Yes",
ifelse(raven_visit_3 !=0, "Yes", "No")))) %>%
mutate(raven_visits = ifelse(raven_visit_3 !=0, 3,
ifelse(raven_visit_2 !=0, 2,
ifelse(raven_visit_1 !=0, 1, 0))))
write.csv(data, "data processed.csv")Bait take by treatment
# Read in processed data
data <- read.csv("data processed.csv") %>%
clean_names() %>%
mutate(date = as.Date(date))
# Plot by treatment type
bait_treat <- ggplot(data, aes(x = as.Date(date),
y = bait_take,
col=treatment,
fill=treatment)) +
geom_smooth() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
ggtitle("All species") +
xlab("Date") +
ylab("Bait take") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(bait_treat)Bait take by treatment and phase
# Read in processed data
data <- read.csv("data processed.csv") %>%
clean_names() %>%
mutate(date = as.Date(date))
# Plot by treatment type
bait_treat_phase <- ggplot(data, aes(x = as.Date(date),
y = bait_take,
col=treatment,
fill=treatment)) +
geom_smooth() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
ggtitle("All species") +
facet_wrap(~phase, scales="free") +
xlab("Date") +
ylab("Bait take") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(bait_treat_phase)# Save the plot as a jpeg
ggsave("figures/bait_treat_phase.jpeg", height=1500, width=2000, unit="px")Percentage bait taken by each species by phase
# Subset to only potential foxes
taken <- data %>%
subset(bait_take != 0) %>%
subset(species != "Not taken") %>%
group_by(species, treatment, phase) %>%
summarise(perc = n() / nrow(data) * 100)
# Plot by treatment type
taken_species_phase_hist <- ggplot(data=taken,
aes(x = species, y = perc,
col=treatment, fill=treatment)) +
geom_col() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
facet_wrap(~phase) +
ggtitle("All species") +
xlab("Species") +
ylab("Percentage of taken baits (%)") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(taken_species_phase_hist)# Save the plot as a jpeg
ggsave("figures/taken_species_phase_hist.jpeg", height=1500, width=2000, unit="px")Plot bait take by rainfall
# Plot by rainfall type
bait_rain <- ggplot(data, aes(x = rainfall_mm,
y = bait_take,
col=treatment,
fill=treatment)) +
geom_smooth() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
ggtitle("All species") +
xlab("Rainfall") +
ylab("Bait take") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(bait_rain)Percentage fox visits by treatment and phase
# Subset to only potential foxes
fox_stats <- data %>%
subset(species == "Fox") %>%
group_by(treatment, phase) %>%
summarise(perc = n() / nrow(data) * 100)
# Plot by treatment type
visits_fox_histogram <- ggplot(data=fox_stats,
aes(y = perc, x = treatment,
col=treatment, fill=treatment)) +
geom_col() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
facet_wrap(~phase) +
ggtitle("Potential fox visits (including unknowns)") +
xlab("Species") +
ylab("Percentage of fox visits (%)") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(visits_fox_histogram)# Save the plot as a jpeg
ggsave("figures/visits_fox_histogram.jpeg", height=1500, width=3000, unit="px")# Subset to only potential foxes
foxes_only <- data %>%
subset(species != "Raven") %>%
subset(species != "Unknown")
# Plot by treatment type
bait_treat_phase <- ggplot(data=foxes_only, aes(x = as.Date(date),
y = bait_take,
col=treatment,
fill=treatment)) +
geom_smooth() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
ggtitle("Foxes vs no take (excluding unknowns)") +
facet_wrap(~phase, scales="free") +
xlab("Date") +
ylab("Bait take") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(bait_treat_phase)Percentage fox visits by treatment and phase
# Subset to only potential foxes
raven_stats <- data %>%
subset(species == "Raven") %>%
group_by(treatment, phase) %>%
summarise(perc = n() / nrow(data) * 100)
# Plot by treatment type
visits_raven_hist <- ggplot(data=raven_stats,
aes(y = perc, x = treatment,
col=treatment, fill=treatment)) +
geom_col() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
facet_wrap(~phase) +
ggtitle("Potential raven visits (including unknowns)") +
xlab("Species") +
ylab("Percentage of raven visits (%)") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(visits_raven_hist)# Save the plot as a jpeg
ggsave("figures/visits_raven_hist.jpeg", height=1500, width=3000, unit="px")# Subset to only potential foxes
ravens_only <- data %>%
subset(species != "Fox") %>%
subset(species != "Unknown")
# Plot by treatment type
bait_raven_treat_phase <- ggplot(data=ravens_only, aes(x = as.Date(date),
y = bait_take,
col=treatment,
fill=treatment)) +
geom_smooth() +
theme_minimal() +
theme(plot.title = element_text(hjust=0.5)) +
ggtitle("Ravens vs no take (excluding unknowns)") +
facet_wrap(~phase, scales="free") +
xlab("Date") +
ylab("Bait take") +
labs(fill="Treatment", col="Treatment")
# Display the plot
print(bait_raven_treat_phase)# Save the plot as a jpeg
ggsave("figures/bait_ravens_treat_phase.jpeg", height=1500, width=3000, unit="px")Bait take by fox visits
# Plot bait take by species visits
bait_visits <- ggplot(data, mapping=aes(x = as.Date(date))) +
geom_smooth(data, mapping=aes(y = as.numeric(bait_take)),
color = "red", fill = "red") +
geom_smooth(data, mapping=aes(y = as.numeric(fox_visits)),
color = "darkorange", fill = "darkorange") +
geom_smooth(data, mapping=aes(y = as.numeric(raven_visits)),
color = "black", fill = "black") +
theme_minimal() +
facet_wrap(~treatment) +
scale_y_continuous(name = "Bait take (%)", limits = c(0, 1),
sec.axis = sec_axis(~., name = "Fox (orange) and raven (black) visits",
breaks = c(0, 1))) +
xlab("Date")
# Display the plot
print(bait_visits)Build generalised linear models for our research questions.
##
## Call:
## glm(formula = bait_take ~ treatment, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.66038 0.02223 29.705 < 2e-16 ***
## treatmentTreatment -0.10124 0.03164 -3.199 0.00142 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.2357498)
##
## Null deviance: 224.02 on 941 degrees of freedom
## Residual deviance: 221.60 on 940 degrees of freedom
## AIC: 1316.1
##
## Number of Fisher Scoring iterations: 2
# Model for bait take by treatment and rainfall
summary(glm(bait_take ~ treatment + rainfall_mm, data=data))##
## Call:
## glm(formula = bait_take ~ treatment + rainfall_mm, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.638871 0.023237 27.493 < 2e-16 ***
## treatmentTreatment -0.100640 0.031505 -3.194 0.00145 **
## rainfall_mm 0.006263 0.002060 3.040 0.00243 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.2337005)
##
## Null deviance: 224.02 on 941 degrees of freedom
## Residual deviance: 219.44 on 939 degrees of freedom
## AIC: 1308.9
##
## Number of Fisher Scoring iterations: 2
##
## Call:
## glm(formula = bait_take ~ treatment + phase, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.69484 0.05679 12.234 < 2e-16 ***
## treatmentTreatment -0.10146 0.03165 -3.205 0.00139 **
## phase -0.01787 0.02710 -0.659 0.50980
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.2358916)
##
## Null deviance: 224.02 on 941 degrees of freedom
## Residual deviance: 221.50 on 939 degrees of freedom
## AIC: 1317.7
##
## Number of Fisher Scoring iterations: 2
##
## Call:
## glm(formula = bait_take ~ treatment + site, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.901039 0.148702 6.059 2.01e-09 ***
## treatmentTreatment -0.102424 0.065102 -1.573 0.116011
## siteF01 0.046330 0.173125 0.268 0.789062
## siteF02 0.201385 0.170912 1.178 0.238991
## siteF03 -0.014532 0.174168 -0.083 0.933523
## siteF04 -0.157455 0.168687 -0.933 0.350858
## siteF05 -0.134257 0.170471 -0.788 0.431161
## siteF06 -0.051039 0.171986 -0.297 0.766718
## siteF07 0.098961 0.170001 0.582 0.560632
## siteF08 -0.003736 0.173070 -0.022 0.982781
## siteF09 0.161306 0.167225 0.965 0.335005
## siteF10 -0.025888 0.171773 -0.151 0.880239
## siteF11 -0.003736 0.173070 -0.022 0.982781
## siteF12 -0.032747 0.169541 -0.193 0.846886
## siteF13 -0.466256 0.169131 -2.757 0.005956 **
## siteF14 -0.560520 0.172711 -3.245 0.001216 **
## siteF15 -0.585677 0.170318 -3.439 0.000611 ***
## siteF16 -0.374846 0.168687 -2.222 0.026524 *
## siteF17 -0.310130 0.170001 -1.824 0.068444 .
## siteF18 -0.480433 0.171773 -2.797 0.005270 **
## siteF19 -0.621452 0.171156 -3.631 0.000298 ***
## siteF20 -0.571342 0.171773 -3.326 0.000916 ***
## siteF21 -0.651039 0.171986 -3.785 0.000164 ***
## siteF22 -0.610447 0.170471 -3.581 0.000361 ***
## siteF23 -0.466256 0.169131 -2.757 0.005956 **
## siteF24 -0.578201 0.169541 -3.410 0.000678 ***
## siteF25 -0.453736 0.173070 -2.622 0.008898 **
## siteF26 -0.434979 0.171773 -2.532 0.011502 *
## siteF27 -0.466256 0.169131 -2.757 0.005956 **
## siteF28 0.201385 0.170912 1.178 0.238991
## siteF29 0.196729 0.171156 1.149 0.250693
## siteF30 0.106147 0.172711 0.615 0.538981
## siteF31 -0.232064 0.172070 -1.349 0.177788
## siteF32 0.058162 0.169541 0.343 0.731635
## siteF34 -0.084329 0.172711 -0.488 0.625480
## siteF35 0.008052 0.170001 0.047 0.962232
## siteF40 -0.666907 0.171156 -3.896 0.000105 ***
## siteF48 -0.751374 0.174168 -4.314 1.78e-05 ***
## siteF52 0.104719 0.167760 0.624 0.532646
## siteF53 0.103414 0.168687 0.613 0.539996
## siteF54 -0.537745 0.170912 -3.146 0.001708 **
## siteF55 -0.091515 0.170949 -0.535 0.592552
## siteF56 0.008601 0.170471 0.050 0.959774
## siteF57 -0.526493 0.167171 -3.149 0.001690 **
## siteF58 -0.810130 0.170001 -4.765 2.20e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1493372)
##
## Null deviance: 223.41 on 937 degrees of freedom
## Residual deviance: 133.36 on 893 degrees of freedom
## (4 observations deleted due to missingness)
## AIC: 924.16
##
## Number of Fisher Scoring iterations: 2
# Do fox visits change with treatment or phase?
summary(glm(fox_visits ~ treatment + phase +
rainfall_mm + maximum_temperature_c, data=data))##
## Call:
## glm(formula = fox_visits ~ treatment + phase + rainfall_mm +
## maximum_temperature_c, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.181054 0.226700 -0.799 0.4247
## treatmentTreatment 0.022190 0.050005 0.444 0.6573
## phase 0.036098 0.053793 0.671 0.5024
## rainfall_mm 0.007035 0.003418 2.058 0.0398 *
## maximum_temperature_c 0.014803 0.009938 1.490 0.1367
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.5886676)
##
## Null deviance: 556.46 on 941 degrees of freedom
## Residual deviance: 551.58 on 937 degrees of freedom
## AIC: 2181.1
##
## Number of Fisher Scoring iterations: 2
# Model for visits against treatment and phase
summary(glm(raven_visits ~ treatment + phase +
rainfall_mm + maximum_temperature_c, data=data))##
## Call:
## glm(formula = raven_visits ~ treatment + phase + rainfall_mm +
## maximum_temperature_c, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.120731 0.101268 1.192 0.233
## treatmentTreatment 0.016744 0.022338 0.750 0.454
## phase -0.002641 0.024030 -0.110 0.913
## rainfall_mm -0.002171 0.001527 -1.422 0.155
## maximum_temperature_c -0.001225 0.004439 -0.276 0.783
##
## (Dispersion parameter for gaussian family taken to be 0.1174665)
##
## Null deviance: 110.37 on 941 degrees of freedom
## Residual deviance: 110.07 on 937 degrees of freedom
## AIC: 662.88
##
## Number of Fisher Scoring iterations: 2
Here we map our responses of bait_take,
fox_visits, and raven_visits across
Ginninderry Conservation Corridor.
# Use Google API to fetch a base map
# ggmap::register_google(key="[add API key here]", write=TRUE)
map <- get_map(location=c(lon=148.9811, lat=-35.2350),
zoom=15, source="google", maptype="satellite", crop=FALSE)# Read in shapefile
gis <- st_read("shapefiles/bait_stations_egg_cta.shp") %>%
fortify() %>%
clean_names() %>%
rename(site=name)## Reading layer `bait_stations_egg_cta' from data source
## `X:\Honours\R2\shapefiles\bait_stations_egg_cta.shp' using driver `ESRI Shapefile'
## Simple feature collection with 43 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 148.9673 ymin: -35.24309 xmax: 148.9927 ymax: -35.22524
## Geodetic CRS: GDA2020
# Map the bait stations
gin_map <- ggmap(map) +
geom_point(data=gis, size=2, aes(x=long, y=lat, col=treatment)) +
theme_minimal() +
xlab("") + ylab("") + labs(col="Treatment")
# Display the plot
print(gin_map)# Combine the bait take df with the GIS df
data <- left_join(data, gis, by="site")
# Generalised linear model
mod <- glm(bait_take ~ lat + long + elevation, data=data)
summary(mod)##
## Call:
## glm(formula = bait_take ~ lat + long + elevation, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.841e+03 7.170e+02 -6.753 2.56e-11 ***
## lat -1.071e+01 3.202e+00 -3.345 0.000857 ***
## long 2.997e+01 4.978e+00 6.020 2.51e-09 ***
## elevation -7.231e-05 6.421e-04 -0.113 0.910369
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1887321)
##
## Null deviance: 222.12 on 930 degrees of freedom
## Residual deviance: 174.95 on 927 degrees of freedom
## (11 observations deleted due to missingness)
## AIC: 1095.7
##
## Number of Fisher Scoring iterations: 2
??annotation_scale
# Map the bait stations
mod_map <- ggmap(map) +
geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
aes(x=long, y=lat, fill=bait_take)) +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),
plot.title=element_text(hjust=0.5),
strip.background=element_rect(fill="#363842"),
strip.text=element_text(colour="white")) +
facet_wrap(~phase) +
scale_fill_viridis_c() +
xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Bait take")
# Display the plot
print(mod_map)# Map the fox visitation
mod_map <- ggmap(map) +
geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
aes(x=long, y=lat, fill=fox_visits)) +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),
plot.title=element_text(hjust=0.5),
strip.background=element_rect(fill="#363842"),
strip.text=element_text(colour="white")) +
facet_wrap(~phase) +
scale_fill_viridis_c() +
xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Fox visits")
# Display the plot
print(mod_map)# Map the fox visitation
mod_map <- ggmap(map) +
geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
aes(x=long, y=lat, fill=raven_visits)) +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),
plot.title=element_text(hjust=0.5),
strip.background=element_rect(fill="#363842"),
strip.text=element_text(colour="white")) +
facet_wrap(~phase) +
scale_fill_viridis_c() +
xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Raven visits")
# Display the plot
print(mod_map)Here we map the bait stations across Ginninderry Conservation Corridor using shapefiles for Ginninderry’s boundaries, paths, roads, and waterbodies.